/*

This code builds data set for model estimation

*/

clear all
set maxvar 10000

*>---------------------- calibrate depreciation -------------------------------<

use "/projects/data/dataSTATA/combined/combined_iv_energy.dta", clear
* Define scalars used in the script in a secondary script for disclosure purposes
do projects/programs/codeSTATA/A_cleaning/5_combine_model_scalars.do

* identify firms that are always surveyed
destring lbdnum, replace
xtset lbdnum year
tsfill

keep if year > 2001 & year < 2008
bys lbdnum: egen count_n = count(firstyear)
keep if count_n == 6
keep year lbdnum ks ck
reshape wide ks ck, i(lbdnum) j(year)
drop if ks2002 == 0
gen n = _n
gen ck2006_ks2007 = ck2006 - ks2007
gen d = .


sum n
* loop over each obs
forvalues i = 1(1)ndepreciation {
preserve

		keep if n == `i'
		sum n
	
		mata: mata clear
		* mata

		mata M = st_data(., .)

		* make sure to preserve order!!
		mata d = polyroots((M[1,15],M[1,11],M[1,9],M[1,7],M[1,5],M[1,2]))

		* mata dReal = select(d,Im(d):==0)
		
		mata dReal = d[1,5]

		mata dReal = Re(dReal)
		
		* mata dReal = dReal[1,1]

		mata st_store(1,16,dReal)

		* end
	keep lbdnum d
	ren d dep
	if `i' > 1 append using `dep'
	tempfile dep
	sa `dep', replace
	
restore

}

merge 1:1 lbdnum using `dep', assert(1 3) keep(1 3) nogen
replace d =  1 - dep
drop dep

sum d


*>--------------------- merge lbd + energy data -------------------------------<

use "/projects/data/dataSTATA/combined/combined_iv_energy.dta", clear

destring lbdnum, replace
xtset lbdnum year

replace ks = . if ks == 0 | ks < 0
gen ks_i = ks 
la var ks_i "Capital Stock, Calculated (total assets ending, 1000 2011 USD)"

** Call scalar to adjust values by

* years before Censuses
bys lbdnum: replace ks_i = (F.ks - ck) / capital_dep if inlist(year,1976, 1981, 1986, 1991, 1996, 2001, 2006, 2011) & !mi(ck)

* years after Censuses
bys lbdnum: replace ks_i = L.ks * capital_dep + L.ck if inlist(year,1978, 1983, 1988, 1993, 1998, 2003, 2008)  & !mi(ck)

* two years before Censuses
bys lbdnum: replace ks_i = (F.ks_i - ck) / capital_dep if inlist(year,1980, 1985, 1990, 1995, 2000, 2005, 2010) & !mi(ck)

* two years after Censuses
bys lbdnum: replace ks_i = L.ks_i * capital_dep + L.ck if inlist(year,1979, 1984, 1989, 1994, 1999, 2004, 2009)  & !mi(ck)

* three years before Censuses
bys lbdnum: replace ks_i = (F.ks_i - ck) / capital_dep if inlist(year,1979, 1984, 1989, 1994, 1999, 2004, 2009) & !mi(ck) & mi(ks_i)

* three years after Censuses
bys lbdnum: replace ks_i = L.ks_i * capital_dep + L.ck if inlist(year,1980, 1985, 1990, 1995, 2000, 2005, 2010)  & !mi(ck) & mi(ks_i)

* four years before Censuses
bys lbdnum: replace ks_i = (F.ks_i - ck) / capital_dep if inlist(year,1978, 1983, 1988, 1993, 1998, 2003, 2008) & !mi(ck) & mi(ks_i)

* four years after Censuses
bys lbdnum: replace ks_i = L.ks_i * capital_dep + L.ck if inlist(year, 1981, 1984, 1991, 1994, 2001, 2004, 2011)  & !mi(ck) & mi(ks_i)

* generate entry/exit flag
gen entrant = 0
replace entrant = 1 if flagb == "Birth"
la var entrant "equals 1 if a plant is an entrant in this year"

gen incumbent = 0
replace incumbent = 1 if flagb == "Continuer"
la var incumbent "equals 1 if a plant neither enters or exits in this year"

gen exit = 0
replace exit = 1 if flagb == "Death"
la var exit "equals 1 if a plant exits in this year"


* interpolate
replace ks_i = . if ks_i < 0 | ks_i == 0
bys lbdnum: ipolate ks_i year, gen(iks_i) epolate 
replace ks_i = iks_i if mi(ks_i) & iks_i > 0
drop iks_i

* do a robustness check of the most excluding these (i.e., drop imputed_k = 1)
gen lks_i = log(ks_i)
reg lks_i c.ltvs##i.bestnaics
predict liks_i, xb
gen iks_i = exp(liks_i)
replace iks_i = ks_i if !mi(ks_i)

la var iks_i "Imputed Capital Stock (total assets ending, 1000 2011 USD)"
gen imputed_k = 0
replace imputed_k = 1 if mi(ks_i) & !mi(iks_i)

gen tot_btu =  btu_int_i_both * tvs * 1000
gen tot_co2 =  co2_int_i_both * tvs * 1000

egen ind_yr_id = concat(bestnaics year)


gen naicscode4 = floor(bestnaics/100)
gen naicscode3 = floor(bestnaics/1000)


save "/projects/data/dataSTATA/combined/combined_model_energy.dta", replace

